#include <math.h>
#include <stdio.h>

/*  ----------- CLASS COMPLEX ----------- */

class Complex {

 public:

  double rl;
  double im;

  Complex() {rl =0; im =0; }
  ~Complex() {}
  Complex(double a, double b) {rl =a; im =b; }
  Complex cplx(double a, double b) {Complex temp; temp.rl =a; temp.im =b; return temp;}
  void operator=(Complex c) { rl = c.rl; im = c.im; }

short operator==(Complex c) {
    if(rl==c.rl && im==c.im) return 1;
    else return 0;
  }
short operator!=(Complex c) {
    if(rl==c.rl && im==c.im) return 0;
    else return 1;
  }
Complex operator+(Complex op) {
    Complex	temp;	
    temp.rl = rl + op.rl;
    temp.im = im + op.im;
    return temp;
}
Complex operator-(Complex op) {
    Complex	temp;	
    temp.rl = rl - op.rl;
    temp.im = im - op.im;
    return temp;
}
Complex operator-(double op) {
    Complex	temp;	
    temp.rl = rl - op;
    temp.im = im;
    return temp;
}
Complex operator*(Complex op) {
    Complex	temp;	
    temp.rl = rl*op.rl - im*op.im;
    temp.im = im*op.rl + rl*op.im;
    return temp;
}
Complex operator*(double op) {
    Complex	temp;	
    temp.rl = rl*op;
    temp.im = im*op;
    return temp;
}
Complex operator/(Complex op) {
    Complex	temp;
    if(op.rl==0 && op.im==0) {
      printf("COMPLEX: Division by zero!\n");
      // exit(1);
    }

    else {
      temp.rl = (rl*op.rl + im*op.im)/(op.rl*op.rl + op.im*op.im);
      temp.im =	(im*op.rl - rl*op.im)/(op.rl*op.rl + op.im*op.im);
    }
    return temp;
}
Complex operator/(double op) {
    Complex	temp;
    temp.rl = rl/op;
    temp.im = im/op;
    return temp;
}
Complex  operator~() {
    Complex	temp;
    temp.rl = rl;
    temp.im = -im;
    return temp;
}

double 	smod() { return rl*rl+im*im;}
double 	mod() { return sqrt(rl*rl+im*im);}
void 	show() { printf(" %le +i %le ",rl,im); }
	
};

inline double cplx_mod_sqr(Complex op) {
  return (op.rl*op.rl+op.im*op.im);
}

inline double cplx_mod(Complex op) {
  return sqrt(op.rl*op.rl+op.im*op.im);
}

inline Complex operator+(double op1,Complex op2) {
  Complex temp;
  temp.rl = op1 + op2.rl;
  temp.im = op2.im;
  return temp;
}

inline Complex operator-(double op1,Complex op2) {
  Complex temp;
  temp.rl = op1 - op2.rl;
  temp.im = - op2.im;
  return temp;
}
inline Complex operator*(double op1,Complex op2) {
  Complex temp;
  temp.rl = op1 * op2.rl;
  temp.im = op1 * op2.im;
  return temp;
}

inline Complex operator/(double op1,Complex op2) {
  Complex temp;
  temp.rl = op1*op2.rl / (op2.rl*op2.rl+op2.im*op2.im);
  temp.im = -op1*op2.im / (op2.rl*op2.rl+op2.im*op2.im);
  return temp;
}

inline Complex cplx_exp(Complex op) {
  Complex temp;
  temp.rl = exp(op.rl) * cos(op.im);
  temp.im = exp(op.rl) * sin(op.im);
			
  return temp;
}

inline Complex	cplx_sqrt(double x) {
  Complex temp;
	
  if(x>=0) {
    temp.rl=sqrt(x);
    temp.im=0.0;
  }
  else {
    temp.rl=0.0;
    temp.im=sqrt(-x);
  }
  return	temp;
}

inline Complex	cplx_sqrt(Complex x) {
  Complex temp;
  double  mod,phase;
	
  if(x.im==0.0)	{
    temp.rl=sqrt(x.rl);
    temp.im=0.0;
  }
  else 	if(x.rl==0.0) {
    if(x.im>0) 	temp.rl=-sqrt(x.im);
    else	temp.rl=sqrt(-x.im);
    temp.im=0.0;
  }
  else {
    mod	= sqrt(sqrt(x.rl*x.rl+x.im*x.im));
    phase = atan(x.im/x.rl);
    temp.rl = mod*cos(phase/2.0);
    temp.im = mod*sin(phase/2.0);
  }
  return temp;
}
